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Abstract 



In various situations in the insurance industry, in finance, in epidemiology, etc., one needs to 
represent the joint evolution of the number of occurrences of an event. In this paper, we present 
a multivariate integer-valued autoregressive (MINAR) model, derive its properties and apply the 
model to earthquake occurrences across various pairs of tectonic plates. The model is an extension 
of (Pedeli and Karlis 2011a) where cross autocorrelation (spatial contagion in a seismic context) is 
considered. We fit various bivariate count models and find that for many contiguous tectonic plates, 
spatial contagion is significant in both directions. Furthermore, ignoring cross autocorrelation can 
underestimate the potential for high numbers of occurrences over the short-term. Our overall 
findings seem to further confirm (Parsons and Velasco 2011). 

Keywords: autoregressive; Granger causality; counts; earthquakes; INAR; multivariate INAR; 
Poisson process 
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1. INTRODUCTION AND MOTIVATION 

1.1 Motivation 

Autoregression in the sense of ARIMA time series models cannot be directly applied to integer 
values for obvious reasons. Thus, most integer-valued autoregressive (INAR) time series models 
are based upon thinning operators such as (Steutel and van Harn 1979) (see also the excellent survey 
of thinning operators by (Wei/3 2008)). Such models have been mainly proposed and investigated 
by (McKenzie 1985) and (Al-Osh and Alzaid 1987) for first order autocorrelation, and by (Du 
and Li 1991) for autocorrelation of order p. (Gauthier and Latour 1994), (Dion, Gauthier and 
Latour 1995) and (Latour 1998) have also investigated a slightly more generalized type of thinning 
operator than (Steutel and van Harn 1979), in models known as generalized INAR (or GINAR). 
The statistical and actuarial literature has multiple successful applications of INAR-type of models 
(see for example (Gourieroux and Jasiak 2004) and (Boucher, Denuit and Guillen 2008) where both 
papers treat car insurance problems). 

In a multivariate setting, the properties of a multivariate INAR (MINAR) model of order 1 
(based upon independent binomial thinning operators) have been derived in (Franke and Subba Rao 
1993) while the multivariate GINAR of order p is presented in (Latour 1997). However, there are 
very few attempts in the literature to estimate and use these types of model^] One notable 
exception is (Pedeli and Karlis 2011a) and (Pedeli and Karlis 20116) who investigated the bivariate 
INAR model of order 1 with Poisson and negative binomial innovations with an application to the 
number of daytime and nighttime accidents. In their papers, the autoregression matrix is diagonal, 
meaning there is no cross-autocorrelation in the counts. 

Insurance policies and earthquake catastrophe (cat) derivatives (such as cat-bonds and cat- 
options) offer protection against earthquake risk in exchange for periodic premiums. Thus, one im- 
portant component in these contracts is the number of earthquakes at various locations. Earthquake 
count models are mostly based upon the Poisson process ((Utsu 1969), (Gardner and Knopoff 1974), 
(Lomnitz 1974), (Kagan and Jackson 1991)), Cox process (self-exciting, cluster or branching pro- 

1 (Heinen and Rengifo 2007) use a Vector Autoregression model for the mean of two Poisson-type of random 
variables. Although the ultimate goal is to represent joint integer-valued random variables, the approach taken is 
very different from multivariate INAR-types of models. 
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cesses, stress-release models (see (Rathbun 2004) for a review), or Hidden Markov Models (HMM) 
(see (Zucchini and MacDonald 2009) and (Orfanogiannaki, Karlis and Papadopoulos 2010)[^J How- 
ever, these models are focused toward a single location whereas seismic risk can also be influenced 
by shocks that occurred at other locations (see e.g. space-time Poisson process in (Ogata 1988), 
(Zhuang, Y. and Vere-Jones 2002) or (Schoenberg 2003)). Thus, one of the purposes of this paper, 
is to propose a bivariate INAR model that accounts for cross-autocorrelation in the counts. From a 
seismological standpoint, that would mean the earthquake count at a given location can be function 
of the past earthquake counts at that site and at another site. These areas can be tectonic plates, 
regions, cities or points on a given geological fault. 

The main objective of this paper is to investigate the effects of seismic space-time contagion (or 
clustering) on various risk management applications using seismological data and a specific model. 
Risk management considerations can be viewed over various time horizons. For example, prices 
of cat-derivatives will be influenced by short-term earthquake risk dynamics because a lack of an 
appropriate earthquake count prediction can mean arbitrage profits or losses may occur on the 
markets. Insurance and reinsurance contracts are managed over a much longer time horizon. 

1.2 Outline of the paper 

(Parsons and Velasco 2011) have confirmed that major earthquakes might have a significant impact 
on the number of earthquakes that occur during the hours following the main shock, but only in 
an area close to the main shock. They do also prove that there is no remote and large earthquakes 
beyond the main shock region. Figure [I] plots the number of quakes following a big one (magnitude 
exceeding 6.5), either within or outside a 2,000 km area from the main shock. One of the aims 
of our paper is to study the dynamics of the number of earthquakes, taking into account spatial 
contagion over tectonic plates. Using plates instead of distance (as in (Parsons and Velasco 2011) 
) allows us to work with multivariate counting processes. 

[Figure 1 about here.] 

The paper is structured as follows. In Section [2j we present the theoretical framework of the 
multivariate INAR of order 1 and the most important results. Some results have already been 
2 For a brief summary of statistical and stochastic models in seismology, see (Vere-Jones 2010). 
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derived in (Franke and Subba Rao 1993) but are presented here for the sake of completeness and 
consistency with our given notation. Additional (theoretical) results in terms of moments, auto- 
covariance functions and predictions are given in this section. In Section [3j we introduce Granger 
causality tests, derived in the given context of BINAR(l) processes, including an interpretation 
of each coefficient in terms of causal effect. Section [4] presents the specific application to the 
bivariate INAR model with cross autocorrelation and Poisson innovations. A Monte Carlo study 
will also illustrate how the maximum likelihood estimators behave (theoretical results are given in 
Section [2j but only in the context of a full cross-correlation matrix) . Finally, Section [5] provides 



various applications of the model with earthquake counts. In subsection 5.3 several BINAR(l) 
processes are fitted over different tectonic plates and magnitudes. We confirm here the conclusions 
of (Parsons and Velasco 2011) claiming that the onset of a large earthquake does not cause other 
large ones at a very long distance. There might be contagion, but it will be between two close 
areas (e.g. contiguous tectonic plates), and over a short period of time (a few hours, perhaps a 



few days, but not much longer). In subsection 5.4, we have also observed that major earthquakes 



will generate several medium-size earthquakes on the same tectonic plate (so called aftershocks). 
Foreshocks were also observed, meaning that medium-size earthquakes might announce the arrival 



of more important earthquakes. To conclude, in subsection 5.5, we compare the sum of counts of 
earthquakes on two plates, assuming that there is - or not - cross correlation between consecutive 
days. 

2. MULTIVARIATE INTEGER- VALUED AUTOREGRESSION OF ORDER 1, MINAR(l) 
As mentioned in (Fokianos 2011), a natural way to define a linear model for counts might be to 
use the Poisson regression to derive an autoregressive process. Let (iVj) denote a count time series, 
and (J-t) the associated filtration. A GARCH-type model can be considered, as in (Ferland, Latour 
and Oraichi 2006) 

p q 

N t \Ft-l ~ V{Xt), where X t = a + s ^a h N t - h + y^f3 k \ t - k - 

h=l k=l 

But one can easily imagine that it could be complicated (and not tractable) to extend such a 
process in higher dimension. An alternative can be to use a thinning operator as in (Al-Osh and 
Alzaid 1987) or (McKenzie 1985). The idea (introduced in (Steutel and van Harn 1979)) is to define 



o as 

p o N = Yi H h Y N if N ^ 0, and otherwise, 

where A" is a random variable with values in N, p G [0,1], and Y"i,l2,-"" are i-i-d- Bernoulli 
variables, independent of N, with P(Yi = 1) = p. Thus po AMs a compound sum of i.i.d. Bernoulli 
variables. Hence, given N, p o N has a binomial distribution with parameters N and p. Based on 
that thinning operator, the integer autoregressive process of order 1 is defined as 

N t =po N t -i +e t =J2 Y i + e *> 

where (et) is a sequence of i.i.d. integer valued random variables. Such process will be called 
INAR(l). Note that such a process can be related to Galton- Watson process with immigration, 
and it is a Markov chain with integer states. As mentioned in (Al-Osh and Alzaid 1987), if (et) 
are Poisson random variables, then (N t ) will also be a sequence of Poisson random variables, and 
the estimation can be done easily using a method of moments estimators or maximum likelihood 
techniques, for p and A = E(e t ). One of the main interest of the thinning operator approach is that 
it can be easily extended in higher dimension, as in (Franke and Subba Rao 1993) (we will also 
provide new results, as well as new interpretations, e.g. in terms of causality). 

2.1 Thinning o operator in dimension d 

As in the univariate case, before defining a multivariate counting process N t := (Aq it , • • • , N^t), we 
need to define a multivariate thinning operator for a random vector N := (Aq, • • • , A^) with values 
in N d . Let P := [pij] be a d x d matrix with entries in [0, 1]. If AT" = (Aq, • • • , A^) is a random 
vector with values in N d , then P o N is a <i-dimensional random vector, with i-th component 

d 

[PoN\ i = Y,p i joX j , 

3=1 

for alH = 1, • • • , d, where all counting variates Y in pij o Xj's are assumed to be independent. 

Note that P o (Q o N) = [PQ] o N. Further, from Lemma 1 in (Franke and Subba Rao 1993), 
E(Po N) = PE(N), and 

E ((P oN)(Po N)') = PE(NN')P' + A, 
with A := diag(VE(iV)) where V is the d x d matrix with entries Pi,j(l — Pi,j)- 
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Definition 2.1 A time series (N t ) with values in N d is called a d-variate INAR(l) process if 

Nt = P o Nt-i + et (1) 

for all t, for some d x d matrix P with entries in [0, 1], and some i.i.d. random vectors et with 
values in N d . 

Remark 2.2 (Pedeli and Karlis 2011a) and (Pedeli and Karlis 2011b) defined bivariate INAR(l) 
processes where matrix P is a diagonal matrix. 

Remark 2.3 (Nt) is a Markov chain with states in N d with transition probabilities 

Tr(n t , nt-!) = F(N t = n t \N t -i = n t -i) (2) 

satisfying 

ir(n t ,n t -i) = ^P(P o = n t - k) ■ P(e = k). 

k=0 

Remark 2.4 Since P has entries in [0,1], using a variant of Perron- Frobenius theorem, there 
exists an eigenvalue K\ of P such that k% > \ni\ for all other eigenvalues of P. And the associated 
eigenvector v\ satistfies V\ > 0. Further, if P has strictly positive entries, pij G (0,1], then 
ki > \m\ and v\ > 0. 

From Remarks |2.3| and \2A\ we can derive sufficient conditions so that there exists a stationary 
MINAR(l) process (based on Theorem 1 and Lemma 2 in (Franke and Subba Rao 1993)). The proof 
is based on the fact that under those assumptions, the Markov chain is irreducible and aperiodic. 
One can prove that is a positive recurrent state, and from Theorem 1.2.2 in (Rosenblatt 1971), 
there exists a strictly stationary solution. 

Proposition 2.5 Let P with entries in (0, 1), such that its largest eigenvalue is less than 1, and 
assume that F(e t = 0) G (0, 1) with Edl^H^) < oo, then there exists a strictly stationary d-variate 
INAR(l) process satisfying Equation Q. 

Lemma 2.6 Let (Nt) denote a stationary d-variate INAR(l) process, with autoregressive matrix 
P with entries in (0, 1), then (Nt) admits a d-variate INMA(oo) representation 

oo 

N t = Y, ph ° e t-h- 

h=0 
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2.2 Maximum likelihood estimation in d-variate INAR(l) processes 

Consider here a finite time series jV = (No, Ni, • • • , N n ), observed from time t = until time 
t = n. The conditional log-likelihood is 

n 

logC(N,e\N ) =J^log?r(JV t _i,JVO (3) 
t=i 

where it is the transition probability of the Markov chain, given by Equation ([2]) . Here parameter 6 
is related to the autoregressive matrix P as well as parameters of the joint distribution of the noise 
process, denoted A. For convenience, assume that A = (Ao, Ai) where Ao are parameters related to 
the innovation process [e-tj margins, and Ai to the dependence among components of the innovation. 
Hence, here = (P, A) on some open sets (0, l)^ 2 x £. From Theorem 2.2 in (Billingsley 1961), 
since (Nt) is a Markov chain, under standard assumptions, we can obtain asymptotic normality of 
parameters. 

Proposition 2.7 Let (Nt) be a d-variate INAR(l) process satisfying stationary conditions, as 
well as technical assumptions (called C1-C6 in (Franke and Subba Rao 1993)), then the conditional 
maximum likelihood estimate of is asymptotically normal, 

Vn(0 -6) 4^(0,E _1 (e)) 1 as n — y oo. 

Further, 

2[log£(iV, 6\N ) - log£QV, 6\N )} A X 2 (d 2 + dim(X)), asn^oo. 

2.3 Autocorrelation matrices for MINAR(l) processes 

Based on the properties obtained in (Franke and Subba Rao 1993) it is possible to derive expressions 
for autocorrelation functions, which is a natural way to describe the dynamics of the process. 

Theorem 2.8 Consider a MINAR(l) process with representation Nt = P ° Nt-i + £t, where 
(et) is the innovation process, with A := E(e^) and A := var(et). Let fj, := E(iV() and ~f(h) := 
cov(N t ,N t _ h ). Then fj, = [I — P]~ X A and for all h € Z, ~f(h) = P h j(0) with 7(0) solution of 
7(0) = P-f(0)P f + (A + A), and where I is the d x d identity matrix. 

Proof 2.9 Since E(Po N) = PE(N), then /x = E(N t ) has to satisfy 

fj, = E(N t ) = E(P o N t -i + e t ) = Pfi + X, i.e. [I - P] ^X. 
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Further 

7 (0) = var(N t ) = E ((P o JVt_i + e t - /x)(P o JV t -i +e t - /x)') 

Since {et) is the innovation process, and from the expression mentioned above (from (Franke and 
Subba Rao 1993)) 

var(N t ) = Pvar(N t -i)P' + A + A 

thus 7(0) satisfies 

7(0) = P7(0)P' + (A + A). 

Finally, 

70) = cov(N t ,N t - h ) = cov((Po N t -i +e t ),N t -h) = cov{{P o (P o N t - 2 + e t -i) + e t ), N t - h ) ■ • • 
eic, so t/iat 

7 (/t) = cow o JV t _ fc + J2 P i ° e t _i, iV t _^ = P h var(N t _ h ) = P h ~/{0), 
since (et) is an innovation process. Q.E.D. 

Remark 2.10 7(0) is a covariance matrix (symmetric) solution of matrix expression 

Z - PZP' = A (= A + A). 

If P was an orthogonal matrix, the term on the left could be related to Lie bracket [Z,P] = 
ZP — PZ. If P was diagonal, we would have obtained expression of (Pedeli and Karlis 2011a). 
Thus, assuming that P can either be orthogonalized or diagonalized will lead to tractable numerical 
algorithm. Another numerical strategy is to seek for a fixed point in equation Z n = PZ n \P' + A 
with some starting value Z$ (e.g. I). This numerical technique will be used in the applications (see 
Section^. 

2.4 Forecasting with MINAR(l) processes 

In order to derive the distribution (or moments) of Nt+h given Nt recall that 

(N t+h , N t ) = (^P h o N t + pl £ t+h-h N t -t^j 
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Proposition 2.11 Let A := E(e t ) and A := var(et), then 

E(N t+h \N t ) = P h N t +(l + p + ... + P h -^ A 

Proof 2.12 The conditional expectation is obtained by recurrence, since the one step ahead condi- 
tional expectation is 

E(N t+1 \N t ) = E(P o N t + e t \N t ) = PN t + A. 

Q.E.D. 

For the conditional variance, it is possible to derive iterative formulas, using recursions. 

Proposition 2.13 Let A := E(e t ) and A := var(e t ), then var(N t+h \N t ) = V h (N t ) where V h (N) 
is defined recursively by Vi(N) = diag{VN) + A and 

V h {N) = E\y h -!(P o N + e)\N] + P^diagiVN) + A](P h ~ 1 y. 

Proof 2.14 The one step ahead conditional variance is 

var(N t+1 \N t ) = var(P o N t + e t |JV t ) = diag(VN t ) + A, 

where V is the dx d matrix with entries Pij(l ~Pi,j)- Then, in order to use a recursive argument, 
we simply have to use the variance decomposition formula, and move one additional step ahead. At 
time t + h, we can write 

var{N t+h \N t ) = E[var{N t+h \ N t+1 )\N t ] + var[E{N t+h \N t+1 )\N t ], 

i.e. 

var(N t+h \N t ) = E[var(N t+h \N t+1 )\N t ] + var[P h - l N t+1 + (l + P + ■ ■ ■ + P h ~ 2 ) X\N t ], 
var(N t+h \N t ) = E[var(N t+h \N t+l )\N t ] + P h - 1 var[N t+1 \N t }(P h - 1 )' . 

Q.E.D. 

In the case where P is diagonal, we obtain as particular case the expressions of (Pedeli and 
Karlis 20116). 
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Corollary 2.15 If P is a diagonal matrix, then 

E(N iit+h \N t ) = p%N itt + [1 +p i4 + ■■■ +p^]X i = p h. Niit + (l^A Xi 

\ 1 Pi,i J 

while 

fl_ v 2h\ (\- V h . l-v 2h \ 

var(N ht+h \N t ) = pj,[l - p^t + r—^f- A* + - — ^ A, 

\ 1 Pi,i J \ 1 Pi,i 1 Pi,i J 

Remark 2.16 If (Nt) is a stationary MINAR(l) process, then, ash — )• oo, K(Nt+h\ Nt) converges 

to E*> ^V = pi--P]-V. 

3. NONDIAGONAL THINNING MATRICES AND GRANGER CAUSALITY 
Based on the concepts introduced in (Franke and Subba Rao 1993), it is possible to get interpre- 
tations of parameters. Based on Granger terminology, N 2 causes Ni at time t if and only if 

E [N lit \Ni,t-i,K2,t-i) * E (Ni, t \N_i,t-i) , 

where = (N 1)0 , ■■■ , Ni, t -i) and N^-i = (N 2}0 , ■■■ , N 2 , t -i). 

Further, N 2 causes instantaneously N\ at time t if 

E (NiMl,t-i,K2,t-l,N 2 , t ) ± E {N ht \Ni,t-i,K2,t-i) ■ 

Thus, as for Gaussian VAR processes, the following interpretation holds (see Section 3.6. in 
(Lutkepohl 2005)) 

Lemma 3.1 Let N t = (Ni jt ,N 2 j) be a bivariate INAR(l) with representation ([I]) 

1. (Ni t) &nd (N%j) are instantaneously related if e is a noncorrelated noise, 

2. (Nit) an d (N 2 ,t) ar e independent, which we denote (Nij) -L (N 2 j), if P is diagonal, i.e. 
pii = p2,i = 0, and eit and e 2y t are independent, 

3. (Ni t t) causes (N 2 ,t) but (N 2 ,t) does not cause (Ni t t), which we denote (Nij) — > (N 2j t), if P 
is a lower triangle matrix, i.e. p 2j i = while pi j2 7^ ; 

4- (N 2 t) causes (Ni t) but (Nit) does not cause (N 2 t), which we denote (Nit) <~ (N 2 t), if P 
is a lower triangle matrix, i.e. pi j2 = while p 2 ,i 7^ 0, 
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5. (JVi,t) causes (Nzj) and conversely, i.e. a feedback effect, which we denote (iVi,t) -H- (N^t), if 
P is a full matrix, i.e. pi 2,P2,l 7^ 



From those characterizations, and Proposition 2.7 it is possible to derive a simple testing 



procedure, based on a likelihood ratio test. For instantaneous causality, we test 

Hq : Ai = A^ against H\ : \\ ^ A^, 

where Ai = \i if and only if margins of the innovation process are independent. 

Corollary 3.2 Let A denote the conditional maximum likelihood estimate of A = (Ao,Ai) in the 
non- constrained MINAR(l) model, and A denote the conditional maximum likelihood estimate of 
\ L = (Ao,A^") in the constrained model (when innovation has independent margins), then under 
suitable conditions, 

2[log£QV,A|iVo) -log£(iV,A |JV )] -4 \ 2 (dim(\) - dim(X ± )), as n -> oo, under H . 
For lagged causality, we test 

H : P eV against Hf.P^V, 

where V is a set of constrained shaped matrix, e.g. V is the set of d x d diagonal matrices for 
lagged independence, or a set of block triangular matrices for lagged causality. 

Corollary 3.3 Let P denote the conditional maximum likelihood estimate of P in the non- constrained 
MINAR(l) model, and P denote the conditional maximum likelihood estimate of P in the con- 
strained model, then under suitable conditions, 

2{log£(N,P\N ) -log£(N,P C \N )} ^> X 2 (d 2 - dim(V)), as n -> oo, under H . 

4. BIVARIATE INAR(l) PROCESS WITH POISSON INNOVATION 
MINAR(d) might appear as tractable models, but the number of parameters can be extremely 
large. In dimension d, the dynamics is characterized by d 2 + dim(A) parameters. The standard 
model for the innovation process would be the multivariate common shock Poisson random vector 
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(see (Mahamunulu 1967) or Section 37.2 in (Johnson, Kotz and Balakrishnan 1997)). Let (ui,I C 
{1, ■ ■ • , d}) be a collection of independent Poisson random variables, and define 

£j := ^ ui 
/C{1,- ,d},i£l 

then e = (ei,-- - ,£d) has a multivariate Poisson distribution. In that case, dim(A) = 2 d . In 
moderate dimension (e.g. d = 10) dim(A) is larger than one thousand, which will not be tractable. 
Thus, for convenience, let us focus on the bivariate INAR(l) process. 

4.1 The bivariate Poisson innovation process 

A classical distribution for e< is the bivariate Poisson distribution, with one common shock, i.e. 

si,t = M 1>t + M , t 
e 2 ,t = M 2 , t + M 0tt 

where M± tt , Mij, and Mq^ are independent Poisson variates, with parameters Ai — if, A2 — <p and 
(p, respectively. In that case, £t = (£i,t,£2,t) has joint probability function 

with Ai,A2 > 0, if G [0, min{Ai, A2}]. See e.g. (Kocherlakota and Kocherlakota 1992) for a 
comprehensive description of that joint distribution. Note that e± t t and £2,* ar e both Poisson 
distributed, with parameter Ai and A2 respectively, and here cov(ei ) t,£2,t) = V 7 - Hence, parameter 
if characterizes independence (or non-independence) of the innovation process. Hence 

and A = 

A 2 



A = 

V 2 / 

and most of the previous expressions can be derived explicitly. 




4.2 BINAR(l) process with Poisson innovation 

For univariate INAR(l) processes, if iVo is assumed to have a Poisson distribution with mean 
X/(l — p), then N t is also Poisson distributed, for all t > 0. But this result does not hold in higher 
dimensions ((Pedeli and Karlis 2011a) noticed that result with diagonal P matrices, and it is still 
true). Nevertheless, it is still possible to derive joint moments of the joint distributions (// and 
7(0)) as well as autocorrelations. 
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Example 4.1 fj, := E(JVj) is given by 

(1 -p2,2)Al +Pl,2A 2 



E(iV 1)t ) = /i! = 

(1 -P2,2) -P2,lPl,2 

(1 -Pl,l)(l ~P2,2) ~P2,lPl,2 



Expressions for 7(0) and 7(1) can be explicitly derived, but from Theorem 2.8 we do have 
matrices based expression that can be used numerically. 

Example 4.2 Auto and cross autocorrelations are given by 

71,2(0) 



corr(N\j, N\ 



2,tj 



v/7i,i (0)72,2(0)' 



{ AT AT \ 71,1 (!) , fAT AT ^ ^i 2 ^ 

corr{Jyij,Jyi,t-i) = tttt ana corr\I\ 2! t, I\ 2t t-i 



corr(N ltt ,N 2t t-i) 



7i,i(0) ' ' ' 72,2(0)' 

7i,2(l) 



V7i,i(0)72, 2 (0)' 



Note that Poisson innovation satisfy technical assumption needed in Proposition 2.7 to insure 
that conditional maximum likelihood estimates converge to a normal distribution as n goes to 
infinity. 

4.3 Maximum likelihood estimation for BINAR(l) with Poisson innovation 

From Equation [3] the conditional likelihood of (P, A, ip) given a sample N = ((Nij, N 2j t), t = 
1,2, ••• ,n) is 

n Ni,t N 2 ,t 

C((P,\,<p);N) = J[ ^2 E *(( N U ~ h, N 2 , t - k 2 ), JV«-i) • P[(ei, t , e 2 ,t) = (fr, fra)] 

t 1 fci fcj A)2 ^2 bivariate Poisson 

with fc x = max{ iV^t - iV ljt _i - N 2 ,t-i, 0} and /c 2 = max{iV 2it - iVi,t-i - N 2 ,t-i, 0}. Here 

7r((ni, n 2 ), JV t _i) = 7Ti(ni, N t -i) • vr 2 (n 2 , iVt-i) 

since given with (et, Nt-i), components of Nt are assumed to be independent (from the definition 
of the multivariate thinning operator o), where m(-,Nt—i) an d ^("i Nt— 1) are convolutions of 
binomial distributions, i.e. for m, n 2 = 0, 1, • • • , N\$-\ + N 2t t-i 

^(ni, jv t _i) = "x: 1 f iVi ' t - i Vbi( i -pi.i) jvi ' t - i - m f iV2 ' t ~ i Vi^a-p^)^- 1 -^-^ 

\ m I ' Vni — ml 
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iVi 



7T 2 (n2,JV t _l) = V f^'WlC 1 -P2,l) Nl - t - 1 ~ m ( W 2 _m (l -p2,2) JV2 ' t - 1 - ( " 2 - m) 

m=0 

Using numerical optimization routines, it is possible to compute (P, A, (p) = argmax{£((P, A, ip); 
and Proposition |2 . 7| insures convergence of that estimator: the conditional maximum likelihood es- 
timates (CMLE) are asymptotically normal and unbiased. 

4.4 Monte Carlo study 

Based on the previous expression of the likelihood, it is possible to run Monte Carlo simulations to 
study the behavior of the estimators on simulated series. This numerical example illustrates with 
two sets of hypothetical parameters how fast is convergence. The two sets of parameters are: 

• pi t % = 0.25, pi : 2 = 0.05, p2 t i = 0.1, p2,2 = 0.4 with Ai = 5, A 2 = 3 and tp = 1; 

• pi t % = 0.25, p\ : 2 = P2,i = 0, p 2 ,2 = 0.4 with Ai = 5, A 2 = 3 and tp = 1. 

The second set of parameters is a special case of the proposed BINAR, which is the diagonal 
BINAR model of (Pedeli and Karlis 20116) and (Pedeli and Karlis 2011a). This will illustrate that 
in some instances such as pi )2 = p 2t i = 0, the CMLE of the multivariate INAR still converges 



to true values (even if Proposition 2.7 insured only convergence in the interior of the support of 
parameters, i.e. (0, l) 4 x (0,oo) 3 , not on borders). 

To perform this experiment, 250 samples of different sizes have been generated. Sample sizes 
of 25, 50, 100, 250, 500, 1000, 5000 and 10000 observations are considered. 

Figure [2] shows the kernel-smoothed^ density function of the distribution of each parameter in 
the first set, over the various sample sizes. Tables [T] and [2] show the mean and standard deviation 
of the parameter values over different sample sizes. One can see that in the first set of parameters, 
the estimates converge quickly to a normal distribution and the bias goes steadily to 0. In the 
second set of parameters, the results are shown in Figure [3] and Tables [3] and [|. O ne sees that even 
though P12 = p 2 ,i = 0, the distribution rapidly concentrates at 0. This indicates that the approach 
is valid. 



3 Although, kernel smoothed densities show some curves in a negative domain, none of estimated parameters was 
negative in the samples. Thus, the negative domain is only due to smoothing. 
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[Figure 2 about here.] 
[Figure 3 about here.] 
[Table 1 about here.] 
[Table 2 about here.] 
[Table 3 about here.] 
[Table 4 about here.] 

5. EMPIRICAL APPLICATION TO EARTHQUAKES 

5.1 Data 

To illustrate the potential of the model for various uses, we apply the proposed BINAR approach 
on earthquake counts of the Earth's tectonic plates. Since the proposed model accounts for se- 
rial correlation and cross-autocorrelation, earthquake counts is an interesting application for the 
following reasons. First, when a mainshock occurs, it provokes many aftershocks, thus creating 
serial correlation. Moreover, the seismic waves travel over a large distance, and may cross different 
tectonic plates, provoking other earthquakes on these other plates (within some time range). This 
is why earthquake counts on contiguous plates should show statistical dependence and cross auto- 
correlation. Given the purpose of the paper, the empirical application is by no means an exhaustive 
seismological analysis of earthquake risk across the planet. From a seismological standpoint, some 
of the results are indicative and further investigation would be required in some aspects of the 
application. 

The data used in the example comes from two sources. First, the limits of each tectonic 
plate come from the Department of Geography of the University of Colorado at Boulder, who 
provide on their website, various shapefiles for use with ArcGlS^ Figure [4] shows the mapping 
of the tectonic plates in the latter reference. The tectonic plates are: North American, Eurasian, 
Okhotsk, Pacific (split in two, East and West), Amur, Indo- Australian, African, Indo-Chinese, 
Arabian, Philippine, Coca, Caribbean, Somali, South American, Nasca and Antarctic. We have 

4 http : //www. Colorado . edu/geography/f oote/maps/ assign/hotspots/hotspots .html 
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decided not to group together the West and East Pacific plates to keep the integrity of the input. 
Secondly, the listings of past earthquakes come from the Advanced National Seismic System (ANSS) 
Composite Earthquake Catalog. Each entry provides the date, time, longitude and latitude, depth 
and magnitude (and its type) of each earthquake. The database spans the time period from 1898 
to 2011, but as mentioned on the ANSS website, many databases have been added between 1898 to 
the mid-1960s. Other factors may have affected the data as well. The addition of seismic stations 
and the technological improvement of seismological instruments may inflate the number of small 
earthquakes in the database. To twart this issue, we focus on earthquakes with a magnitude of at 
least 5 (M > 5), and we used a subset of the data. To find the most appropriate cutoff date in the 
data, we used statistical tests of changes in structures (F test (Chow test), from (Andrews 1993) 
or (Zeileis, Kleiber, Kramer and Hornik 2003) for implementation issues). Based on these tests, 
events between January 1st, 1965 up to March 30th, 2011 have been kept in the sample. The 
total number of earthquakes is approximately 70 000. Finally, for M > 6 earthquakes, the time 
period considered is January 1st, 1992 up to March 30th, 2011, which amount to approximately 
3000 events. 

The proposed BINAR model will be largely used to investigate first-order autocorrelation and 
cross-autocorrelation in earthquake counts, which is equivalent to measuring the degree of the first 
order type of space and time contagion. To do this, earthquake counts have been computed at 
several frequencies. Time ranges of 3, 6, 12, 24, 36 and 48 hours have been considered to count the 
number of earthquakes. 

[Figure 4 about here.] 

5.2 Quality of fit 

The proposed BINAR model encompasses various models as well. When, pij = 0,i,j = 1,2 and 
if = 0, there is no serial correlation, no cross autocorrelation, and both series are independent. 
Those are two independent Poisson noises. When ip / 0, the Poisson noises are dependent. When 
Pi, 2 = P2,i = and ip = 0, there is serial correlation but both series are independent. Those are 
equivalent to two univariate INAR processes. Finally, when = p2,i = and ip / 0, we find 
the diagonal BINAR model of (Pedeli and Karlis 2011a). In the latter model, there is no cross 
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autocorrelation. Thus, in this section we compare the fit of those four models, along with the 
proposed BINAR approach, on each of the 136 possible pairs of tectonic plates. 

Table [5] shows the results of the likelihood ratio test (LRT) of 2 dependent Poisson noises and 
of 2 independent INARs, both compared with 2 independent Poisson noises. Each column shows 
descriptive statistics for various sampling frequencies. The meaning of each row is as follows: mean 
LRT across the 136 combinations of tectonic plates, standard deviation of LRTs, along with various 
quantiles (50%, 75%, 90%, 95%, 97.5%) and the proportion of combinations that is statistically 
significant. Thus, the value of the LRT provides an idea of how useful the added feature really is. 
The upper part of Table [5] shows that dependence in the noises is important for about 10% of the 
combinations of plates (most of them are contiguous). However, serial correlation is a much more 
important feature, even though noises are independent. This should have been expected given 
earthquake mechanics. Thus, the LRT shows that autocorrelation is one important component 
when analyzing earthquake counts at these sampling frequencies. 

Sampling frequency also influences the degree of autocorrelation and cross autocorrelation. 
When the sampling frequency is h hours, all earthquakes on two tectonic plates will count towards 
cp in the time interval [0, h] hours. All earthquakes occurring in the time interval [h, 2h] will help 
find first degree autocorrelation and cross autocorrelation i.e. it should appear in the P matrix. 
Finally, all earthquakes occurring after 2h hours, would be accounted for if a second or third degree 
BINAR was considered. Thus, when h increases, p should become more significant. This is what 
we observe in the upper panel of Table [5] Even though the percentage of combinations that are 
significant very slightly increase, the mean LRT and higher percentiles of LRT tend to grow more 
importantly. Note that only a few combinations are such that tectonic plates are contiguous. 

Table [6] shows similar computations for the diagonal BINAR model over the independent INAR 
model (contribution of p), and for the proposed BINAR model over the diagonal BINAR model 
(contribution of p% 2 and P2,i)- One sees that 6-13% of combinations of tectonic plates show a 
significant dependence in the noise (upper part of Table [6]) . Those are in large part contiguous 
plates, which explains the rather low percentage. For other combinations, plates are too far apart 
and independent INAR models are often sufficient. When we compare the proposed BINAR model 
to the diagonal model (lower part of Table pi) , we find that 6-9% of the combination of plates show 
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cross autocorrelation. In other words, the non-diagonal terms in the P matrix, i.e. and P2,i, 
are both statistically different from zero. In most of the cases, the combination of plates that had 
a significant fit to both models, are contiguous. We further investigate some of those in the next 
subsection. 

5.3 Analysis of pairs of tectonic plates 

In this section, we take a look at specific pairs of tectonic plates to observe parameters and interpret 
them. We will look at four different combinations of plates, which are all closely related to Japan. 
The Japanese area is one of the most seismically active regions of the world, being at the limit 
of 4 tectonic plates. Table [7] shows the (CMLE) parameter estimates at four different sampling 
frequencies, for the four combinations of plates. The 8th line of each panel displays the LRT over 
the diagonal BINAR model. A value of more than 5.99 is significant at a level of 95%, meaning 
that both cross autocorrelation terms are significant. The last two lines show the unconditional 
mean number of earthquakes per period on each plate. 

The first and second order moments estimators for the Okhotsk and West Pacific plates are 
given in Table [8] 

Let us focus on the Okhotsk and West Pacific plates, where the results are shown at the bottom 
left part of Table [7j and assume the sampling frequency is 24 hours. Thus, the daily number of 
earthquakes on the Okhotsk plate is explained by three sources: the number of earthquakes on the 
previous day on both plates, and a random noise effect. When no earthquake was observed on both 
plates on a given day, a Poisson r.v. with mean 0.16 earthquake per day will generate seismicity 
on the next day. The probability of observing one or more earthquakes by noise only is 15% in this 
case. 

The interest of the proposed BINAR model lies in the representation of the spatial contagion 
effect between tectonic plates. Suppose that n earthquakes were observed on the Okhotsk plate on 
a given day, while m earthquakes were observed on the West Pacific plate on that same day. The 
number of earthquakes on the Okhotsk plate the next day will be the result of the convolution of 
a binomial(n, 0.0817) (autocorrelation of order 1), a binomial(m, 0.028) (cross autocorrelation of 
order 1) and a Poisson noise (mean 0.162). Thus, on average, the number of earthquakes on the 
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next day on the Okhotsk plate will be 

0.0817ra + 0.028m + 0.162. 
Under a diagonal model estimated with CMLE (but not provided in Table [7]), that quantity is 

0.0922n + 0.1748 

which ignores the contribution of the West Pacific plate's earthquakes. With n > and m > 1, 
the diagonal model will understate the potential number of earthquakes on the Okhotsk plate. 
When to is relatively large, which is one of the cases we are interested in risk management, that 
understatement can be important. For example, if 3 earthquakes are observed on the West Pacific 
plate, and only one on the Okhotsk plate, then we have an average of 0.3277 earthquake on the 
latter plate with the proposed BINAR model, compared with 0.267 with the diagonal model. 

Further, suppose a case where n = and to > 1, in a context where we want to compare the 
mean number of earthquakes on the Okhotsk plate using the diagonal and proposed models. In 
that situation, the mean number of earthquakes in the proposed model is roughly 16% larger]^] than 
the mean number of earthquakes in the diagonal model, for each additional earthquake we observe 
on the West Pacific plate (i.e. 16% times to). Picking other sets of plates, even if the LRT is much 
smaller and still significant, will lead to similar analyses (Okhotsk and Amur at 12- or 24-hour 
frequency is one example among others). 

One may be tempted to directly compare values of and p2,i and conclude that earthquakes on 
one plate provokes earthquakes on the other. However, the gross values of pi 5 2 and p2,i ar e obviously 
influenced by the number of earthquakes on each plate. For example, at the 24-hour frequency, 
pi 2 = 0.028 and p2,i = 0.106 so that one may mistakenly pretend that Okhotsk earthquakes 
generally provoke earthquakes on the West Pacific plate, and not the converse. However, there 
are approximately 3 times more earthquakes on the West Pacific plate than on the Okhotsk plate, 
meaning that pi_2 has to be lower to compensate for the higher counts on the second plate. Thus, if 
one is interested in determining if earthquake counts on one plate determine the other, one should 
perform Granger causality tests. 

5 02 8 7+° 8 162 - 1 is approximately ^§ = 0.16m. 
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[Table 5 about here.] 
[Table 6 about here.] 
[Table 7 about here.] 
[Table 8 about here.] 

5.4 Foreshocks and aftershocks 

As another application of the model, we illustrate the relationship between medium-size earthquakes 
(i.e. 5 < M < 6) and large earthquakes (M > 6). Using the proposed BINAR model in this context 
will help understand how the size of a set of earthquakes at a given time period can help predict 
the size of future earthquakes. Most of the time, large earthquakes (mainshocks) are followed 
by aftershocks, which are usually smaller (medium-size or smaller). The inverse, in which case a 
medium-size earthquake may announce a larger earthquake, is usually less likely, but still regularly 
observed. Figure [5] illustrates this relationship between foreshocks, mainshocks and aftershocks. 

[Figure 5 about here.] 

As a first exercise, we have fitted the same five models, that is the proposed BINAR model, the 
diagonal BINAR model, independent INARs, dependent Poisson noises, and independent Poisson 
noises. According to the LRT, the fit of the diagonal BINAR model over independent Poisson 
noises is statistically significant for all tectonic plates, at all sampling frequencies. This is also the 
case when the diagonal BINAR is compared to independent INAR models. Finally, for all but a 
few tectonic plates and/or sampling frequencies, the diagonal BINAR model has a very significant 
fit over dependent Poisson noises. 

Thus, for this application, we would like to measure if cross autocorrelation is important, i.e. 
if earthquake size on a given period helps explain future earthquake sizes. Table [9] shows the LRT 
for the proposed BINAR model over the diagonal model, for various sampling frequencies. A value 
larger than 5.99 means that ^ and p2,i 7^ 0, implying that large earthquakes are followed by 
medium-size earthquakes, and the opposite also holds. This is indeed the case in the large majority 
of tectonic plates, although this relationship clearly gets weaker when the sampling frequency goes 



21 



from 3- hours to 48-hours (last row of the table). This should have been expected given Omori's 
law, which explains the temporal decay of aftershock rates. 



Table 10 shows the CMLE parameter estimates for all plates at the 12-hour frequency. The 
last two columns provide the unconditional mean number of medium (5 < M < 6) and large 
(M > 6) earthquakes. With a 12-hour sampling frequency, only the Coca and Somali plates 
have an unsignificant LRT at a level of 95%, meaning that for 15 plates, cross autocorrelation 
is important. Thus, one should not directly compare values of p\^ and p% i since only Granger 
causality tests will provide the true significance of contagion between the two sets of data. 

Let us illustrate the impact of cross autocorrelation for a given tectonic plate. Assume that 
on the Okhotsk plate, which seats beneath part of Japan, there is a large earthquake in the prior 
12-hour period (and no medium-size earthquake). Then, cross autocorrelation will be the most 
important component of the mean number of earthquakes in the next period. Indeed, the expected 
number of medium earthquakes in the next period is 0.2444 + 0.0780 = 0.3224 and cross autocorre- 
lation will account for more than the two thirds of the total expectation. One can compare the size 
of pi t 2 and P2.1 with the noise components (As) and observe that the ratio is much larger in this 
section than in Section [5 .3| Thus, cross sectional effects are a key element in this context. Finally, 
the ratio of the expected number of M > 6 earthquakes over 5 < M < 6 earthquakes is on average 
(across plates) approximately 10 which is consistent with Gutenberg and Richter's law. 

[Table 9 about here.] 

[Table 10 about here.] 

5.5 Risk management 

In many risk management applications, such as the computation of premiums and reserves or the 
pricing of catastrophe derivatives, the total loss amount over a given area, region, or city is what 
matters most. One important driver of the total loss amount, is the total number of earthquakes over 
the area in question for various time horizons [0, T] . In this section, we compare the distribution of 
the sum of the number of earthquakes over a given area, for the diagonal and the proposed BINAR 



models. We do so for pairs of tectonic plates (see Section 5.3) where the LRT was statistically 
significant, otherwise the two models are too similar. 
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Two sets of tectonic plates are analyzed: (1) Okhotsk and West Pacific plates (Japan is at the 
limits of the West Pacific, Okhotsk, Philippine and Amur plates) and (2) South American and Nasca 
plates (which holds the South American continent and the West Coast of South America (Chile for 
example) is located at the limit of these two plates). Two extreme scenarios are generated. In the 
first set of plate, we assume that 23 earthquakes were observed on the Okhotsk tectonic plate and 
46 were observed on the West Pacific plate (this is indeed what happened in the last 12 hours of 
March 10th, 2011). In the second set of plates, we assume that 24 earthquakes were observed on 
the South American plate, whereas only 3 were observed on the Nasca plate (this is what occurred 
on the second half of February 27th, 2010). Using 100 000 paths of a bivariate diagonal INAR 
and the proposed bivariate INAR models, we have computed the total number of earthquakes that 
occurred on both plates (of a given set), on the next T days (T = 1,3,7, 14 and 30). The results 



are shown in Table 11 The left (right) panel focuses on the first (second) set of tectonic plates. 

J-q ) for various values of n. 



The numbers shown are P ^$3fc=i C^l,fc + ^2,k) > n 

One sees that the diagonal model really understates the number of earthquakes in the following 
days, especially in the tails. For example, in the first set of plates (Okhotsk and West Pacific), the 
probability of having a total of at least 20 earthquakes in the next day is 6.7% with the proposed 
model, whereas it is 0.7% with the diagonal model; it is a ten-fold increase. This increase is all due 
to the non-diagonal terms in the P matrix as it accounts for the cross auto-correlation between 
earthquake counts. A less dramatic increase is observed in the second set of plates (South America 
and Nasca). For example, the probability of having a total of at least 7 earthquakes over a week on 
both plates is 39.6% in the diagonal model whereas this probability is 44% in the proposed model. 
As expected, over the long-term, both processes converge to their equilibrium and the effect of the 
initial conditions seem to dissipate. 



We now suppose that with both sets of plates, no earthquake occurred on a given day. Table 12 



shows the results of P [Ylk=i {^i,k + ^2,fc) > 71 -^oj f° r T = 14 and 30 days. For smaller T values, 
the probabilities generated by the two models are very similar since it takes a lot of time to develop 
earthquakes and thus to observe cross-sectional effects. For the given T values, the probabilities 
are very similar for both models, with a slightly fatter tail for the proposed model in the first set 
of tectonic plates. In the second set of plates, the probabilities are too close to be able to conclude 
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of any difference. 

In summary, we have also run different scenarios on other sets of plates and it confirms that the 
effect of the non-diagonal terms in the P matrix is to generate fatter tails in the sum of the number 
of earthquakes. This is very useful for short-term risk management applications such as the pricing 
of earthquake bonds and other derivatives. An underestimation of the number of earthquakes could 
mean arbitrage opportunities if the market model has a similar behavior to the diagonal model. 

[Table 11 about here.] 

[Table 12 about here.] 

6. CONCLUSION 

In this paper, we confirm the conclusion of (Parsons and Velasco 2011) claiming that very large 
earthquakes do not necessarily cause large ones at a very long distance. There might be contagion, 
but it will be within two close areas (e.g. contiguous tectonic plates), and over a short period of 
time (a few hours, perhaps a few days, but not much longer). Nevertheless, not taking into account 
possible spatial contagion between consecutive periods may lead to large underestimation of overall 
counts. In the context of foreshocks, mainshocks and aftershocks, we have also observed that major 
earthquakes might generate several medium-size earthquakes on the same tectonic plate, and also 
foreshocks might announce possible large earthquakes. 
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List of Figures 

1 Number of earthquakes (magnitude exceeding 2.0) per 15 seconds, following a large 
earthquake (of magnitude 6.5) either close to the main shock (less than 2,000 km) or 
far away (more than 2,000 km). Counts were normalized so that the expected number 
of earthquakes before is 100 in the two regions. Plain lines are spline regressions, 
either before or after the main shock) [29] 

2 Distribution of estimators pi t i, pi,2, P2,l, Pi,ii -Vi.) ^2 and <p, as a function of the 
sample size n, case of non-diagonal P matrix [30] 

3 Distribution of estimators pi,i, Pi,2, P2,i, P2,2, -Vlj ^2 and tp, as a function of the 
sample size n, case of diagonal P matrix [31] 

4 The 17 tectonic plates (North American, Eurasian, Okhotsk, Pacific (split in two, 
East and West), Amur, Indo- Australian, African, Indo-Chinese, Arabian, Philippine, 
Coca, Caribbean, Somali, South American, Nasca and Antarctic) [32] 

5 Number of earthquakes (magnitude exceeding 2.0) per 15 seconds, following a large 
earthquake (of magnitude 6.5), normalized so that the expected number of earth- 
quakes before is 100. Plain lines are spline regressions, either before or after the 
main shock) [33] 
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Figure 1: Number of earthquakes (magnitude exceeding 2.0) per 15 seconds, following a large 
earthquake (of magnitude 6.5) either close to the main shock (less than 2,000 km) or far away 
(more than 2,000 km). Counts were normalized so that the expected number of earthquakes before 
is 100 in the two regions. Plain lines are spline regressions, either before or after the main shock). 
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Figure 2: Distribution of estimators pi t i, Pi,2, P2,i, P2,2, Ai, A2 and <p, as a function of the sample 
size n, case of non-diagonal P matrix. 



30 



p[1 .1]: True value is 0,25 




p[2,1 ]: True value is 



|j[2,2]. True value is 0.4 
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Figure 3: Distribution of estimators pip, P2,i, P2,2 ; Ai, A2 and (p, as a function of the sample 
size n, case of diagonal P matrix. 
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Figure 4: The 17 tectonic plates (North American, Eurasian, Okhotsk, Pacific (split in two, East 
and West), Amur, Indo- Australian, African, Indo-Chinese, Arabian, Philippine, Coca, Caribbean, 
Somali, South American, Nasca and Antarctic). 
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Figure 5: Number of earthquakes (magnitude exceeding 2.0) per 15 seconds, following a large 
earthquake (of magnitude 6.5), normalized so that the expected number of earthquakes before is 
100. Plain lines are spline regressions, either before or after the main shock). 
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Sample size n 


Pi,i 


Pl,2 


P2,l 


P2,2 


Ai 


A 2 


<P 


25 


21.27% 


12.15% 


14.54% 


34.24% 


4.8023 


2.9854 


1.1057 


50 


23.18% 


9.97% 


11.73% 


38.62% 


4.8538 


2.9850 


0.9704 


100 


23.29% 


6.63% 


10.11% 


39.93% 


5.0075 


3.0081 


1.0034 


250 


23.95% 


5.52% 


10.47% 


39.56% 


5.0531 


2.9801 


0.9774 


500 


24.57% 


6.19% 


10.18% 


39.49% 


4.9704 


3.0318 


1.0162 


1000 


24.93% 


5.02% 


10.09% 


39.84% 


5.0044 


3.0040 


0.9843 


5000 


24.88% 


4.96% 


9.95% 


39.92% 


5.0097 


3.0086 


0.9952 


10000 


25.06% 


4.98% 


10.08% 


39.94% 


4.9969 


2.9972 


1.0036 


True value 


25% 


5% 


10% 


40% 


5 


3 


1 



Table 1: Mean parameter values - First parameter set 
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Sample size n 


Pi,i 


Pl,2 


P2,l 


P2,2 


Ai 


A 2 


<P 


25 


16.87% 


15.43% 


14.77% 


19.42% 


1.2599 


1.1232 


0.9243 


50 


13.28% 


11.24% 


10.58% 


11.98% 


1.0497 


0.9038 


0.7480 


100 


9.25% 


7.89% 


7.99% 


8.35% 


0.7336 


0.6218 


0.5344 


250 


6.04% 


5.24% 


5.52% 


5.33% 


0.5043 


0.4063 


0.3701 


500 


3.94% 


4.44% 


3.82% 


3.72% 


0.3601 


0.3106 


0.2516 


1000 


2.94% 


3.22% 


2.74% 


2.55% 


0.2587 


0.2144 


0.1813 


5000 


1.37% 


1.54% 


1.17% 


1.19% 


0.1148 


0.1002 


0.0766 


10000 


0.92% 


1.00% 


0.83% 


0.84% 


0.0841 


0.0660 


0.0568 



Table 2: Standard deviation of parameter values - First parameter set 
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Sample size n 


Pi,i 


Pl,2 


P2,l 


P2,2 


Ai 


A 2 


<P 


25 


20.87% 


11.18% 


8.15% 


36.00% 


4.6948 


2.7182 


1.1375 


50 


23.09% 


6.10% 


5.29% 


38.16% 


4.8105 


2.7352 


1.0503 


100 


23.22% 


5.47% 


3.20% 


39.66% 


4.8545 


2.7997 


1.0095 


250 


25.08% 


2.59% 


2.09% 


39.24% 


4.8667 


2.9008 


1.0442 


500 


24.76% 


1.98% 


1.38% 


39.91% 


4.9131 


2.9150 


1.0333 


1000 


24.93% 


1.42% 


1.00% 


40.22% 


4.9382 


2.9211 


0.9906 


5000 


24.84% 


0.72% 


0.38% 


40.00% 


4.9740 


2.9743 


1.0017 


10000 


25.05% 


0.44% 


0.34% 


39.92% 


4.9738 


2.9820 


1.0018 


True value 


25% 


0% 


0% 


40% 


5 


3 


1 



Table 3: Mean parameter values - Second parameter set 
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Sample size n 


Pi,i 


Pl,2 


P2,l 


P2,2 


Ai 


A 2 


<P 


25 


17.48% 


16.80% 


11.65% 


17.46% 


1.2515 


0.9555 


0.8950 


50 


13.42% 


11.42% 


7.75% 


12.09% 


0.9161 


0.7196 


0.6621 


100 


10.05% 


7.83% 


4.83% 


8.07% 


0.7286 


0.4765 


0.4751 


250 


5.71% 


4.19% 


3.08% 


5.37% 


0.4263 


0.3293 


0.3125 


500 


4.07% 


2.73% 


2.04% 


3.62% 


0.3101 


0.2191 


0.2048 


1000 


2.82% 


2.00% 


1.36% 


2.48% 


0.1981 


0.1605 


0.1624 


5000 


1.33% 


0.97% 


0.61% 


1.14% 


0.0992 


0.0702 


0.0696 


10000 


0.88% 


0.68% 


0.53% 


0.78% 


0.0635 


0.0519 


0.0446 



Table 4: Standard deviation parameter values - Second parameter set 
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Likelihood ratio test - Dependent Poisson over independent Poisson 





3 hours 


6 hours 


12 hours 


24 hours 


36 hours 


48 hours 


Mean 


6.9755 


8.9167 


9.7735 


9.9896 


10.6302 


10.4286 


Stdev 


67.8029 


76.2997 


83.3619 


85.4670 


91.0044 


88.9381 


50% 


0.0000 


0.0000 


0.0000 


0.0000 


0.0000 


0.0087 


75% 


0.0315 


0.4181 


0.5951 


1.0178 


1.0533 


1.2933 


90% 


1.6329 


3.6041 


4.3060 


3.7918 


4.1652 


4.7956 


95% 


5.1601 


6.6625 


11.2303 


11.7174 


12.8974 


10.5802 


97.5% 


9.9646 


16.9966 


18.8452 


22.7211 


24.8961 


19.5306 


% > 3.84 


7.35% 


10.29% 


12.50% 


10.29% 


11.03% 


12.50% 


Likelihood ratio test - independent INARs over 


independent Poisson 




3 hours 


6 hours 


12 hours 


24 hours 


36 hours 


48 hours 


Mean 


1215.16 


1150.81 


1036.95 


864.04 


557.15 


542.15 


Stdev 


1084.90 


1082.05 


964.69 


798.20 


483.23 


456.70 


50% 


851.01 


781.32 


735.33 


561.29 


399.63 


416.24 


75% 


1630.64 


1497.45 


1415.05 


1173.66 


819.89 


745.15 


90% 


2979.68 


3030.71 


2678.53 


2147.04 


1423.80 


1319.45 


95% 


3227.93 


3170.56 


2837.07 


2308.66 


1515.82 


1435.19 


97.5% 


3551.20 


3589.28 


3213.06 


2580.33 


1761.71 


1650.87 


% > 5.99 


100.00% 


100.00% 


100.00% 


100.00% 


100.00% 


100.00% 



Table 5: Likelihood ratio test (1) independent Poisson vector (with A ^ 0) over independent Poisson 
variables (2) two independent INAR processes versus two independent Poisson variables 
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Likelihood ratio test - diagonal BINAR over independent INARs 





3 hours 


6 hours 


12 hours 


24 hours 


36 hours 


48 hours 


Mean 


3.3744 


4.8765 


5.4843 


6.3690 


9.1950 


8.0845 


Stdev 


26.9264 


36.7864 


42.7106 


52.6144 


72.6791 


63.8500 


50% 


0.0048 


0.0391 


0.0955 


0.0171 


0.2013 


0.0337 


75% 


0.5303 


0.5109 


0.7107 


0.8735 


1.1663 


1.0702 


90% 


1.8276 


2.7837 


4.4992 


5.0022 


4.2161 


4.3279 


95% 


4.9423 


4.7359 


8.3814 


9.9875 


10.7675 


8.4055 


97.5% 


9.2657 


15.0495 


13.3836 


13.4121 


24.5204 


16.5784 


% > 3.84 


6.62% 


8.82% 


12.50% 


13.24% 


10.83% 


10.83% 


Likelihood ratio test - proposed BINAR over 


diagonal BINAR 




3 hours 


6 hours 


12 hours 


24 hours 


36 hours 


48 hours 


Mean 


4.9409 


4.1814 


3.7077 


4.0927 


2.8335 


3.5242 


Stdev 


30.4265 


25.5655 


23.8860 


22.6459 


12.6545 


15.1581 


50% 


0.9720 


0.5379 


0.3533 


0.4504 


0.3631 


0.5931 


75% 


2.6904 


2.3514 


1.6824 


2.4026 


2.0943 


2.3357 


90% 


5.2349 


5.2194 


4.0567 


4.7533 


4.4141 


5.3640 


95% 


9.9654 


7.8503 


6.4292 


8.5268 


6.5271 


8.7984 


97.5% 


15.7839 


15.6885 


12.0481 


11.9986 


11.9423 


18.2802 


% > 5.99 


8.09% 


8.82% 


7.35% 


8.82% 


6.67% 


9.02% 



Table 6: Likelihood ratio test (1) diagonal BINAR over two independent INAR processes (2) 
proposed BINAR over diagonal BINAR, with Poisson innovation. 
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I ) 1 J 

relates 


Okhotsk (#1) vs 


Philippine (#2) 


Okhotsk (#1) 


vs. Amur 


(#2) 


Params /Frequency 


3 hours 


12 hours 


24 hours 


48 hours 


3 hours 


12 hours 


24 hours 


48 hours 


Pi,i 


7.44% 


9.45% 


10.38% 


12.81% 


7.44% 


9.44% 


10.30% 


12.75% 


Pi, 2 


0.61% 


0.60% 


1.15% 


0.00% 


0.35% 


0.83% 


3.06% 


2.31% 


P2,l 


0.00% 


0.00% 


0.00% 


0.00% 


0.16% 


0.42% 


0.44% 


0.40% 


P2,2 


3.87% 


5.83% 


8.52% 


8.80% 


4.68% 


6.44% 


8.67% 


10.59% 


\ 

M 


0.0222 


0.0868 


0.1711 


0.3358 


0.0223 


0.0871 


0.1720 


0.3348 


A 2 


0.0156 


0.0612 


0.1187 


0.2368 


0.0032 


0.0122 


0.0237 


0.0466 


<? 


0.0000 


0.0001 


0.0000 


0.0021 


0.0000 


0.0003 


0.0009 


0.0024 


Lnl (over uiag.j 


4.1106 


1.6737 


3.6011 


0.0000 


2.8874 


9.4113 


9.5405 


4.0631 


Uncond. mean (#1) 


0.0241 


0.0963 


0.1926 


0.3852 


0.0241 


0.0963 


0.1927 


0.3852 


Uncond. mean (#2) 


0.0162 


0.0650 


0.1298 


0.2596 


0.0034 


0.0134 


0.0269 


0.0538 


Plates 


Okhotsk (#1) vs. 


West Pacific (#2) 


Okhotsk (#1) vs. 


IndoChinese (#2) 


Params /Frequency 


3 hours 


12 hours 


24 hours 


48 hours 


3 hours 


12 hours 


24 hours 


48 hours 


Pi,i 


6.12% 


7.18% 


8.17% 


10.13% 


7.45% 


9.46% 


10.36% 


12.83% 


Pi, 2 


1.85% 


2.85% 


2.80% 


3.13% 


0.02% 


0.28% 


0.24% 


0.10% 


P2,l 


5.84% 


7.56% 


10.60% 


9.74% 


0.22% 


0.40% 


0.00% 


0.75% 


P2,2 


10.71% 


13.52% 


15.52% 


15.67% 


6.71% 


10.29% 


11.58% 


13.68% 


Al 


0.0214 


0.0818 


0.1620 


0.3132 


0.0223 


0.0863 


0.1710 


0.3344 


A 2 


0.0576 


0.2212 


0.4261 


0.8539 


0.0767 


0.2948 


0.5818 


1.1326 


<P 


0.0012 


0.0098 


0.0269 


0.0739 


0.0002 


0.0015 


0.0046 


0.0099 


LRT (over diag.) 


352.5998 


275.2342 


257.0215 


157.0995 


0.2839 


3.0150 


1.2208 


0.6136 


Uncond. mean (#1) 


0.0241 


0.0963 


0.1926 


0.3852 


0.0241 


0.0963 


0.1926 


0.3852 


Uncond. mean (#2) 


0.0661 


0.2643 


0.5285 


1.0570 


0.0823 


0.3290 


0.6580 


1.3155 



Table 7: Estimation of parameters for counts of earthquakes on several tectonic plates, Okhotsk vs. 
Philippine; Okhotsk vs. Amur; Okhotsk vs. West Pacific; and Okhotsk vs. Indo-Chinese plates. 
Includes a likelihood ratio test (null: P is a diagonal matrix). 
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Plates 


Okhotsk (#1) vs. 


West Pacific (#2) 


Params / Frequency 


3 hours 


12 hours 


24 hours 


48 hours 


E(JVi, t ) 


0.024 


0.096 


0.192 


0.385 


E(iV 2it ) 


0.065 


0.264 


0.528 


1.057 


var(iV lit ) 


0.022 


0.084 


0.167 


0.326 


var(iV 2 ,t) 


0.060 


0.239 


0.466 


0.934 


coT(N ltt ,N 2 ,t) 


0.038 


0.079 


0.110 


0.150 


cor(iV M ,iVi )t _i) 


0.062 


0.075 


0.086 


0.109 


cor(jV 2)t ,iV2, t _i) 


0.108 


0.138 


0.162 


0.165 


cor(7Vi )t ,JV 2)t -i) 


0.033 


0.053 


0.055 


0.068 



Table 8: First and second order moments, fi and 7(0) and cross-lagged correlations p(0) and p(l), 
for counts on two plates, Okhotsk vs. West Pacific. 
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Plate name 


3 hours 


6 hours 


12 hours 


24 hours 


36 hours 


48 hours 


North American 


60 


9470 


31 


3554 


21 


6005 


16 


5618 


15 


9287 


6 


0150 


Eurasian 


3 


4172 


2 


4860 


17 


5732 


1 


3990 


8 


4201 


1 


0743 


Okhotsk 


135 


5948 


109 


1666 


109 


3060 


113 


5049 


36 


3703 


52 


2677 


East Pacific 


37 


2827 


50 


2991 


32 


2566 


19 


9613 


19 


1339 


4 


4437 


West Pacific 


101 


3846 


96 


3865 


110 


2205 


109 


0303 


62 


4744 


81 


3029 


Amur 


12 


9162 


17 


2652 


7 


7396 


4 


0498 


10 


2767 


15 


0012 


Indo-Australian 


303 


2257 


233 


9429 


169 


1037 


124 


9187 


75 


0355 


48 


7183 


African 


35 


0197 


11 


1661 


15 


9194 


12 


7146 


28 


1233 


9 


0930 


Indo-Chinese 


63 


2515 


29 


9391 


49 


5970 


64 


0781 


29 


8289 


45 


1555 


Arabian 


4 


5921 


4 


5763 


12 


3768 


3 


1358 


2 


1765 





1744 


Philippine 


9 


2969 


21 


3144 


20 


1805 


15 


5858 


18 


9310 


17 


5329 


Coca 


12 


6070 


15 


1147 


3 


1709 


8 


2198 


5 


3469 


9 


0246 


Caribbean 


20 


4764 


24 


5509 


21 


2256 


3 


4367 


7 


6112 


2 


4771 


Somali 





2432 


5 


1726 


3 


2162 





0039 





1625 





0392 


South American 


81 


8145 


58 


0135 


50 


4781 


55 


8060 


77 


3867 


58 


4621 


Nasca 


76 


8393 


38 


6514 


20 


2549 


17 


6382 


8 


6659 


11 


1903 


Antarctic 


2 


9275 


9 


2410 


9 


1584 


4 


7911 


1 


8339 





9290 


Average LRT 


56 


5786 


44 


6260 


39 


6105 


33 


8139 


23 


9827 


21 


3471 



Table 9: Likelihood Ratio Test for the proposed BINAR model over the diagonal model, for vari- 
ous sampling frequencies, when ATi jt denotes the number of medium size earthquakes (magnitude 
between 5 and 6) during period t, and Nij denotes the number of large earthquakes (magnitude 
exceeding 6) during period t. 
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Plate name 




Pl,l 




Pl,2 




P2,l 




P2,2 




Ai 




A 2 




<p 


Pl,2/P2,l 


Mean<6 


Mean > 6 


North American 





05633 





11372 





00444 





01027 





0844 





0091 





0028 


25.63 





0906 





0096 


Eurasian 





01348 





14431 





01082 





00006 





0181 





0011 





0002 


13.34 





0185 





0013 


Okhotsk 





11224 





24445 





00995 





01951 





0780 





0104 





0033 


24.57 





0910 





0115 


East Pacific 





07950 





12959 





00385 





00025 





2631 





0285 





0075 


33.64 





2900 





0296 


West Pacific 





15688 





21797 





00642 





01163 





1995 





0212 





0084 


33.93 





2426 





0231 


Amur 





00931 





09470 





00676 





02041 





0107 





0024 





0008 


14.02 





0110 





0025 


Indo- Australian 





19749 





24562 





01079 





03095 





4039 





0490 





0225 


22.76 





5205 





0564 


African 





03906 





13683 





00204 





00716 





0564 





0054 





0014 


67.20 





0595 





0056 


Indo-Chinese 





09744 





16198 





00563 





00956 





2501 





0236 





0080 


28.76 





2816 





0254 


Arabian 





04026 





24457 





00347 





00009 





0167 





0007 





0001 


70.51 





0176 





0008 


Philippine 





03630 





09681 





00864 





03512 





0536 





0055 





0012 


11.20 





0563 





0062 


Coca 





06228 





04115 





00155 





00534 





0439 





0069 





0020 


26.61 





0471 





0070 


Caribbean 





03080 





26310 





00001 





00009 





0083 





0008 





0004 


31969 





0088 





0008 


Somali 





02325 





02809 





00384 





00000 





0284 





0012 





0001 


7.32 





0291 





0013 


South American 





13661 





12043 





01141 





01507 





1384 





0160 





0046 


10.55 





1628 





0181 


Nasca 





11426 





13361 





00307 





01442 





0378 





0034 





0013 


43.49 





0433 





0036 


Antarctic 





02875 





03897 





00879 





00153 





0548 





0056 





0010 


4.43 





0567 





0061 



Table 10: CMLE estimators for the proposed BINAR model, for 12-hour frequency, when N±j 
denotes the number of medium size earthquakes (magnitude between 5 and 6) during period t, and 
N\j denotes the number of large earthquakes (magnitude exceeding 6) during period t. 
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Diagonal model (Okhotsk and West Pacific) 


Diagonal model (South American and Nasca) 


n / days 


1 day 


3 days 


7 days 


14 days 


30 days 


n / days 


1 day 


3 days 


7 days 


14 days 


30 days 


5 


0.9680 


0.9869 


0.9978 


0.9999 


1.0000 


2 


0.8489 


0.9166 


0.9757 


0.9981 


1.0000 


10 


0.5650 


0.7207 


0.8972 


0.9884 


0.9999 


5 


0.2708 


0.4321 


0.6965 


0.9277 


0.9988 


15 


0.1027 


0.2270 


0.4978 


0.8548 


0.9985 


7 


0.0685 


0.1628 


0.3959 


0.7655 


0.9906 


20 


0.0067 


0.0277 


0.1308 


0.4997 


0.9752 


10 


0.0035 


0.0192 


0.1041 


0.4108 


0.9334 


25 


0.0003 


0.0018 


0.0170 


0.1684 


0.8588 


15 


0.0000 


0.0002 


0.0033 


0.0547 


0.5885 


30 


0.0000 


0.0001 


0.0014 


0.0319 


0.5965 


20 


0.0000 


0.0000 


0.0000 


0.0031 


0.1873 


40 


0.0000 


0.0000 


0.0000 


0.0002 


0.1034 


25 


0.0000 


0.0000 


0.0000 


0.0001 


0.0290 


50 


0.0000 


0.0000 


0.0000 


0.0000 


0.0041 


30 


0.0000 


0.0000 


0.0000 


0.0000 


0.0030 


Proposed model (Okhotsk and West Pacific) 


Proposed model (South American and Nasca) 


n / days 


1 day 


3 days 


7 days 


14 days 


30 days 


n / days 


1 day 


3 days 


7 days 


14 days 


30 days 


5 


0.9946 


0.9977 


0.9997 


1.0000 


1.0000 


2 


0.8780 


0.9321 


0.9805 


0.9979 


1.0000 


10 


0.8344 


0.9064 


0.9712 


0.9970 


1.0000 


5 


0.3323 


0.4888 


0.7331 


0.9362 


0.9990 


15 


0.3638 


0.5288 


0.7548 


0.9479 


0.9995 


7 


0.0990 


0.2034 


0.4410 


0.7913 


0.9921 


20 


0.0671 


0.1573 


0.3616 


0.7256 


0.9917 


10 


0.0082 


0.0309 


0.1271 


0.4435 


0.9386 


25 


0.0053 


0.0246 


0.0970 


0.3815 


0.9357 


15 


0.0000 


0.0004 


0.0056 


0.0688 


0.6145 


30 


0.0002 


0.0023 


0.0151 


0.1268 


0.7646 


20 


0.0000 


0.0000 


0.0001 


0.0039 


0.2099 


40 


0.0000 


0.0000 


0.0001 


0.0038 


0.2335 


25 


0.0000 


0.0000 


0.0000 


0.0001 


0.0380 


50 


0.0000 


0.0000 


0.0000 


0.0001 


0.0221 


30 


0.0000 


0.0000 


0.0000 


0.0000 


0.0036 



Table 11: Empirical evolution of P (El=i (-^1,* + N 2,k) > n Fo) for various values of n (per line) 
and T (per column), on two plates (Okhotsk vs. West Pacific and South American vs. Nasca), 
either for a diagonal P matrix - on top - or for a full matrix - below. 
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Diagonal model 


(Okhotsk and West Pacific plates) 


(South American and Nasca plates) 


n / days 


14 days 


30 days 


n / days 


14 days 


30 days 


5 


0.9495 


1.0000 


3 


0.9096 


0.9991 


10 


0.5281 


0.9943 


5 


0.6710 


0.9904 


15 


0.1138 


0.9125 


7 


0.3733 


0.9506 


20 


0.0121 


0.6302 


10 


0.0962 


0.7770 


25 


0.0008 


0.2791 


12 


0.0296 


0.5883 


30 


0.0000 


0.0774 


15 


0.0038 


0.2996 


35 


0.0000 


0.0137 


20 


0.0001 


0.0498 


40 


0.0000 


0.0018 


25 


0.0000 


0.0036 


45 


0.0000 


0.0002 


30 


0.0000 


0.0002 


50 


0.0000 


0.0000 


35 


0.0000 


0.0000 


Proposed model 


(Okhotsk and West Pacific plates) 


(South A 


merican 


and Nasca plates) 


n / days 


14 days 


30 days 


n / days 


14 days 


30 days 


5 


0.9444 


0.9999 


3 


0.9061 


0.9990 


10 


0.5211 


0.9927 


5 


0.6662 


0.9899 


15 


0.1261 


0.9033 


7 


0.3754 


0.9497 


20 


0.0139 


0.6242 


10 


0.0990 


0.7745 


25 


0.0007 


0.2877 


12 


0.0317 


0.5856 


30 


0.0000 


0.0870 


15 


0.0049 


0.2986 


35 


0.0000 


0.0177 


20 


0.0001 


0.0513 


40 


0.0000 


0.0024 


25 


0.0000 


0.0044 


45 


0.0000 


0.0003 


30 


0.0000 


0.0002 


50 


0.0000 


0.0001 


35 


0.0000 


0.0000 



J-q ) for various values of n and two 



Table 12: Empirical evolution of P (J2k=i ( N hk + N 2,k) > n 
time horizon T, on two plates (Okhotsk vs. West Pacific and South American vs. Nasca), either 
for a diagonal P matrix, or for a full matrix. 
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